#ifndef _HYBROID_H_
#define _HYBROID_H_
#define MaxNumNeighbors 4


typedef struct node {
  double uu, vv;		/* parametric coords */

  double mass;			/* nodal mass */
  double damp;			/* nodal damping */
  vector refer;			/* reference component wrt model frame*/
  vector displ;			/* displacement component wrt model frame */
  vector ddispl;                /* time derivative of displacement */
  vector posit;			/* position component wrt to global frame */
  vector force;			/* nodal applied force */
  vector rforce;                /* applied force in reference orientation */
  vector veloc;			/* nodal velocity */
  
  double sinu;			/* precomputed sin(u) */
  double cosu;			/* precomputed cos(u) */
  double sinv;			/* precomputed sin(v) */
  double cosv;			/* precomputed cos(v) */

  int node_id;			/* Zhennan: id for node */

  int num_neighbors;		/* number of neighbors */
  struct node *neighbors[MaxNumNeighbors]; /* indices of neighboring nodes */
  struct node **many_neighbors; /* special case if needed:
					 indices of many neighbors */

  struct elem **elements;    /* array of pointers to corresponding elements */
  int max_elem_num;         /* number of corresponding elements */
  int *corr_number;        /* array of corresponding number of node in 
                             **elements array */
  vector average_normal;

} node;


typedef struct link {
  node *p1, *p2;		/* pointers to linked nodes */
} link;

typedef struct elem {
  node *vertices[MaxNumNeighbors]; /* array of element nodes */
  int nomnodes;                  /* number of elements */
  int elemtype;                  /* 0 for rectangular elms */
                                 /* 1 for triangular for south pole */
                                 /* 2 for triangular for north pole */
                                 /* 3 for triangular in the mesh */
} elem;

typedef struct range_datapoint {
  vector position;		/* position of datapoint in space */
  double u, v;			/* point of influence of data force */
  node *nearest_node;		/* pointer to nearest node on model */
  double strength_force;	/* strength of data force */
  } range_datapoint;


typedef struct hybroid {
  int usize, vsize;		/* parametric dimensions */
  double ubegin, uend, ustep;   /* ranges and spacings */
  double vbegin, vend, vstep;   /* of uniform parametric grid */
  int numnodes, numlinks;	/* number of nodes and links */
  int numelems;                 /* number of elements */
  node *nodes;			/* array of nodes */
  elem *elems;                  /* array of elements */
  double s400, s401, s402,      /* components of symmetric stiffness matrix */
         s403, s411, s412,      /* for 4 node element */
         s413, s422, s423,
         s433;
  double s300_1, s301_1, s302_1,/* components of symmetric stiffness matrix */
         s311_1, s312_1, s322_1;/* for 3 node element for south pole*/
  double s300_2, s301_2, s302_2,/* components of symmetric stiffness matrix */
         s311_2, s312_2, s322_2;/* for 3 node element for north pole*/
  double s300_3, s301_3, s302_3,/* components of symmetric stiffness matrix */
         s311_3, s312_3, s322_3;/* for 3 node element for middle element*/

  double n4_1[3], n4_2[3], n4_3[3], n4_4[3]; /* matrices of shape functions */
                     /* calculated at the 4 Gauss points for 4 node element */
  double n31_1[2], n31_2[2], n31_3[2],
         n32_1[2], n32_2[2], n32_3[2]; /* matrices of shape functions */
        /* calculated at the 3 Radau points for the two types of 3 node elms */
  
  link *links;			/* array of links */
  
  double w0;			/* deformation magnitude resistance */
  double w1;			/* deformation gradient resistance */
  
  double sc, sc_pr, a1, a2, a3, e1, e2; /* superquadric parameters */
  double dsc, da1, da2, da3, de1, de2; /* derivatives of above parameters */
  double tapx, tapy, bend1, bend2, bend3; /* tapering and bending params */
  double dtapx, dtapy, dbend1, dbend2, dbend3; /* deroivatives of above para */
  
  double mass;			/* global mass */
  double damp;			/* global damping */
  double dt;			/* time step */
  double sdt;                   /* time step for parameter changes */

  int numconstr;               /* number of constraints, 4 is maximum */
  struct constraint *constr[4]; /* pointers to the four constraints maximum */
  int suqid[4];               /* denotes suq number in constraint,i.e.1 or 2 */
  int flag;                    /* flag stating the visit of hybroid */
  int fflag;           /* denotes the application of a force on a constraint */

  vector translation;		/* global translation */
  vector dtranslation;          /* time derivative of global translation */
  quaternion orientation;       /* total orientation */
  quaternion dorientation;      /* time derivative of orientation */ //Mingchen, B in thesis
  double rotation[3][3];	/* global rotation */
  vector net_applied_force, net_applied_torque;	/* total forces on model */
  vector net_damping_force, net_damping_torque;

  range_datapoint *range_data;	/* pointer to range dataset */
  int numdatapoints;		/* number of points in dataset */
  
  int *image;			/* image potential */
  int image_xsize, image_ysize;	/* image sizes */
  double strength_image_force;
}hybroid;


typedef struct velocity {
  struct constraint *constr; /* constraint which causes */
                             /* these velocities to occur */
  vector velf0;
  vector velf1;
  vector velf2;
  vector velf3;
}velocity;

typedef struct constraint {
  int number;      /* global number of the constraint  */

  hybroid *h1;
  int point1;
  velocity vel1_pr[1][4]; /* assume four constraints maximum per object */
                         /* velocity of this constraint point under loads */
                 /* on itself and on the other constraints of this hybroid */
  velocity vel1[1][4]; /* assume four constraints maximum per object */
                            /* velocity of this constraint point under loads */
                   /* on itself and on the other constraints of this hybroid */


  hybroid *h2;
  int point2;
  velocity vel2_pr[1][4]; /* assume four constraints maximum per object */
                         /* velocity of this constraint point under loads */
                   /* on itself and on the other constraints of this hybroid */
  velocity vel2[1][4]; /* assume four constraints maximum per object */
                         /* velocity of this constraint point under loads */
                  /* on itself and on the other constraints of this hybroid */

  vector fc; /* force acting on this constraint on hybroid h1*/
  vector b;
}constraint;



#endif _HYBROID_H_